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I FOREWORD 

This  report  was  prepared  by  E.  W.  Turner,  Aerospace  Engineer  in 
; the  Loads  and  Response  Prediction  Group  of  the  Structural  Mechanics 

[ ( Division  of  the  Air  Force  Flight  Dynamics  Laboratory  at  Wright- 

Patterson  Air  Force  Base,  Ohio.  The  work  described  herein  is  a part 
of  the  Air  Force  Systems  Command  exploratory  development  program  to 
predict  aircraft  response  to  atmospheric  turbulence.  The  work  was 
directed  under  Project  1367,  "Structural  Integrity  for  Military 
Aerospace  Vehicles,"  Task  136701,  "Structural  Flight  Loads  Data." 

f j 

This  report  covers  work  done  in  the  period  from  October  1975  to 
August  1976.  This  manuscript  was  released  by  the  author  in  September 

[ 1976  for  publication  as  an  AFFDL  Technical  Report. 
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SECTION  I 
INTRODUCTION 

The  traditional  Power  Spectral  Density  (PSD)  analysis  for 
calculating  the  response  of  an  aircraft  to  continuous  turbulence  has  evolved 
through  a continual  effort  to  account  more  precisely  for  the  unpredictable 
nature  of  turbulence.  Although  few  will  question  that  the  PSD  analysis 
does  a considerably  better  job  of  representing  the  true  nature  of 
turbulence  than  the  discrete  gust  model,  the  traditional  PSD  analysis 
has  been  questioned  because  it  under  predicts  the  presence  of  turbulence 
in  the  higher  intensities. 

This  report  is  an  effort  to  view  the  traditional  PSD  analysis  and 
the  reported  deficiences  in  their  proper  perspective.  Hence  the  entire 
approach,  depicted  in  Figure  1,  and  the  basic  assumptions  are  reviewed. 

The  deficiencies  are  then  discussed,  and  finally,  the  most  recent  efforts 
to  characterize  turbulence  are  presented. 
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SECTION  II 

ATMOSPHERIC  TURBULENCE  MODELS 

One  day  after  a particularly  rough  flight  in  gusty  winds,  Wilbur 
supposedly  remarked  to  Orville,  "Once  we  get  above  the  trees,  we  won't 
have  to  worry  about  this  bumpiness."  Similarly,  the  prevailing  thought 
until  not  too  long  ago  was  that  flights  above  the  clouds  would  be  in 
still  air.  The  discovery  of  the  jet  stream,  temperature  inversions, 
mountain  waves,  and  other  types  of  turbulence  proved  these  theories  to 
be  false. 

Above  the  region  in  which  the  atmosphere  behaves  as  a boundary 
layer,  the  numerous  types  of  turbulence  are  frequently  divided  into  two 
categories:  (1)  convective  turbulence  in  and  around  clouds,  particularly 
thunderstorms,  and  (2)  CAT,  which  is  an  acronym  for  Clear  Air  Turbulence. 
Below  the  cloud  bases,  surface  heating  may  result  in  direct  buoyant 
convection  and  many  small  eddies  that  are  perceived  as  mild  CAT.  Above 
a field  of  cumulus  clouds,  a regular  "chop"  may  persist  where  there  Is 
significant  gradients  in  the  wind  velocity  with  altitude.  Some  of  the 
more  violent  types  of  CAT  are  associated  with  mountain  and  lee  waves 
which  may  have  rotors  (vortices)  embedded  in  them  at  regular  intervals. 
Flights  in  the  vicinity  of  the  jet  stream  and  near  the  tropopause  can  be 
particularly  turbulent,  depending  on  meteorological  conditions.  Some  of 
the  most  violent  turbulence  of  all  exists  in  thunderstorms  and  squall 
lines,  particularly  in  areas  of  heavy  precipitation.  The  mechanism 
(Reference  1)  for  generating  this  turbulence  is  convection  which  produces 
rising  and  falling  columns  of  air  ringed  by  toroids  of  extreme  vorticity. 

The  mechanism  of  turbulence  is  such  a varied  and  complicated  process 
that  statistics  offer  the  only  manageable  method  to  handle  the  gust 
design  problem. 

1 . DISCRETE  GUSTS 

For  the  purpose  of  design  for  gust  encounter,  a description  (or 
model)  of  atmospheric  turbulence  is  required.  For  simplicity,  an  analytical 
model  is  desirable.  Early  models  characterized  turbulence  as  a discrete 
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gust  and  the  "sharp  edged  gust"  is  of  this  type.  Neglecting  lift  growth, 
the  maximum  response  of  a rigid  aircraft  to  a sharp  edged  gust  is: 


(1) 


where: 


nQ  = load  factor  in  steady  level  flight  =1.0 

n = load  factor  in  turbulence 

12  2 
q = dynamic  pressure  = 2P0ve»  lbs/ft 

S = wing  area,  ft2 

= lift  curve  slope 

W = gross  weight,  lbs 

Vg  = airspeed  (equivalent),  ft/sec 

U = gust  velocity  (equivalent),  ft/sec 

e 2 4 

p = air  density  at  sea  level,  lb  sec  /ft 


This  model  evolved  ^nto  a "one-minus-cosine"  gust  in  an  effort  to  use  a 
more  realistic  waveform.  The  factor  1^  was  added  to  account  for  this 
waveform  and  also  to  account  for  the  finite  length  of  time  required  for 
lift  to  develop  the  following  gust  encounter,  which  is  referred  to  as 
"lift  growth." 


n « n©  + 


Pa  vt  Cla  ud«  "w 
2W/S 


(2) 


By  substituting  Vkeas  = 1.689  Vg  and  pQ  = .002377,  this  equation  becomes 
the  so-called  "gust  loads  formula"  which  has  been  used  in  design  extensively. 


where: 


498 

Ude 

v 

keas 

W/S 


n 


v«*»  cla  Ud« 
"°  498  W/S 


2.0 

(1.689)  (0.002377) 

derived  gust  velocity,  ft/sec 

knots  equivalent  airspeed 

2 

wing  loading,  lbs/ft 


(3) 

(4) 
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Values  for  the  somewhat  arbitrary  factor  1^  are  given  in  terms  of 
aircraft  mass  ratio  jj  for  aircraft  at  subsonic  airspeeds  as: 

K.  . -MBS-  ( 

■ 5.3+p. 

and  for  aircraft  at  supersonic  airspeeds  as: 


" 6.95  + /*' 

where 

u s 2-W/A  .. 
^ 9 * * 


and  o 

g = 32.2  ft/sec^ 

c = mean  aerodynamic  chord,  feet 
p = air  density  at  flight  altitude,  lb  sec2/ft4 


Design  values  of  gust  velocity  for  use  in  the  gust  loads  formula 
were  established  following  analysis  of  many  hours  of  VGH  (velocity- 
load factor-altitude)  data.  Essentially,  the  gust  loads  formula  was 
solved  in  reverse  for  measured  occurrences  of  load  factor  peaks,  to 
determine  what  gust  velocity  from  a 25-chord  gust  would  produce  the 
observed  peak.  Accepted  design  values  for  the  derived  gust  velocity 
are  given  in  Figure  2.  These  values  are  specified  in  paragraph 
25.341  of  the  Federal  Aviation  Regulations. 


25  38 


25  50  66 

Ude  GUST  VELOCITY  (ft./»ec.) 

Figure  2.  Derived  Gust  Velocity  for  Gust  Loads  Formula 
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where 

VL  = limit  speed  for  the  basic  high  drag  (landing)  configuration 
= maximum  continuous  speed  in  level  flight 
Vq  = gust  penetration  speed 


i 

! 

! 


! 

! 

I 


I 


L 


2.  C0NTiN:"'JS  TURBULENCE  MODELS 

The  25-chord  wavelength  for  the  one-minus-cosine  gust  was  selected 
because  this  wavelength  historically  couples  with  the  vertical  trans- 
lation and  short-period  pitch  modes  of  the  rigid  airplane  to  produce 
the  largest  induced  load  factor.  With  the  development  of  higher  speed 
aircraft  with  increased  flexibility,  the  possibility  of  other  wavelengths 
coupling  with  flexible  modes  to  oroduce  significant  response  is  Increased. 
In  an  effort  to  consider  all  of  the  possible  wavelengths,  and  the  magnitude 
of  the  gust  associated  with  those  wavelengths,  investigators  turned  to 
power  special  techniques  for  aircraft  design.  This  analysis  is  based 
on  a model  which  characterizes  turbulence  as  a stationary  random  process; 
considerable  time  and  effort  was  expended  to  collect  sufficient  data  to 
establish  this  model. 

As  part  of  an  Air  Force  sponsored  "Critical  Air  Turbulence  (CAT)" 
investigation,  instruments  were  installed  on  fixed  towers  to  record 
true  gust  velocities  at  low  altitudes,  and  instrumented  aircraft  were 
flown  to  collect  flight  data  under  various  meteorological  conditions 
and  altitudes.  The  data  which  must  be  collected  in  flight  in  order 
to  calculate  time  histories  of  gust  velocity  consist  of  continuous  re- 
cording of  airflow  direction,  airspeed,  and  aircraft  motions.  Three 
axis  airflow  directions  were  recorded  from  flow  vanes  mounted  on  a 
boom  on  the  nose  of  the  test  aircraft.  The  measured  parameters  which 
define  the  aircraft  motion  consist  of  three  axis  accelerations  (ax, 
ay’  az}  an<*  two  axis  rotat1on  rates  (pitch  and  yaw  rates).  These 
parameters  are  used  in  conjunction  with  the  equations  of  motion  of 
the  aircraft  to  convert  angle  of  attack  of  the  boom  mounted  vanes 
to  gust  velocities  and  to  remove  the  contributions  from  the  aircraft 
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motion.  The  following  equation  illustrates  this  for  the  case  of  vertical 
gust  velocity: 

Wfl  * (av  - 6 ) VT  ♦ JL%  6 + fat  dt  (8) 

where 

Wg  = true  vertical  gust  velocity 

ay  = angle  of  attack  of  boom  mounted  vane 

0 = pitch  angle  of  attack 

6 = pitch  rate  of  aircraft 
Vj  = true  velocity  of  aircraft 

lx  = distance  from  vane  to  center  of  gravity  of  aircraft 
az  = vertical  acceleration  of  aircraft 

All  of  the  parameters  are  measured  from  the  mean  values  taken  over 
the  complete  record. 

Although  more  precise  equations  are  sometimes  used,  it  is  con- 
venient to  assume  that  the  aircraft  and  boom  are  completely  rigid, 
and  that  the  angles  involved  in  the  calculations  of  gust  velocities  are 
small  enough  that  the  accuracy  is  not  lost  by  assuming  small -angle 
theory.  Variations  in  upwash  at  the  vane  location  on  the  boom  are 
either  neglected,  or  are  accounted  for  in  the  calibration  of  the 
vanes.  Some  sample  time  histories  of  vertical  and  lateral  gust  velocity 
appear  in  Figure  3 (Reference  2). 

3.  POWER  SPECTRA 

The  Power  Spectral  Density  (PSD)  of  a time  varying  function  X(t) 
is  a real  function  which  shows  how  the  mean  square  value  of  X(t)  is 
distributed  with  frequency.  The  ordinate  of  the  PSD  curve  at  a par- 
ticular frequency,  u,  is  the  mean  square  value  of  that  part  of  X(t) 
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Figure  3.  Sample  Time  Histories  of  Vertical  Gust  Velocity 
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whose  frequency  is  within  an  infinitely  narrow  band  centered  about 
the  frequency  w.  This  may  be  expressed  as: 

$(«)  = i‘m,t  -~r—  / (9) 

Aui — ► o Tflw 

T — ► oo 


where 

<5( ti)  = power  spectral  density  function  of  X(t)  in  units  of 
(ft/sec) (rad/sec) 


T = duration  of  X(t)  in  seconds 


X(t,(jMoj)  = component  of  X(t)  within  frequency  band  u + Au> 

T 


The  total  area  under  the  PSD  curve  is  the  mean  square  value  of  X(t) 
taken  over  all  frequencies. 

In  order  to  remove  the  effects  of  variations  in  flight  speed 
from  one  analysis  to  another,  power  spectra  of  the  atmosphere  are 
calculated  on  the  basis  of  a spacial  frequency  defined  as: 

n = -f-  do) 

VT 

where 

ft  = spacial  frequency  in  radians  per  foot 

w = circular  frequency  as  observed  by  the  aircraft  in  radians 
per  second 

VT  = true  velocity  of  the  aircraft  in  feet  per  second 


The  relation  between  the  original  and  the  transformed  power  spectra 
is  given  by: 


$(«)  = 


01) 
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Although  it  is  possible  to  calculate  the  power  spectra  directly  from 
the  definition,  it  is  more  convenient  to  first  calculate  the  auto- 
correlation function,  R(t),  for  the  complete  record,  and  then  obtain 
the  power  spectrum,  d(n),  by  performing  a Fourier  transform  of  the 
autocorrelation  function  using  the  Cooley-Tukey  Fast  Fourier  Transform 
(FFT)  technique. 

R(r)  = /x( t)x(t  + r)dt  (12) 

0 

and 

Ha)’~k  JR(T>‘‘aTdT  (13> 

—CD 


Typical  power  spectra  of  three  different  meteorological  conditions 
are  showr.  in  Figure  4 (Reference  2).  The  height  of  the  curve  is  a 
measure  of  the  intensity  of  the  turbulence  at  particular  frequencies, 
and  these  show  a characteristic  decrease  in  intensity  with  increase  in 
frequency.  The  square  root  of  the  area  under  the  power  spectra  curve  is 
a measure  of  the  overall  rms  gust  velocity,  a.  Since  the  power  spectra, 
in  theory,  extends  from  zero  frequency  to  infinite  frequency,  and  in 
practice  there  is  some  maximum  frequency  which  is  considered  called 
the  cutoff  frequency,  uic,  there  are  two  possible  values  of  o.  In 
practice,  only  the  area  under  the  measured  power  spectrum  is  included 
in  the  rms  value,  even  though  the  total  rms  value,  if  all  of  the  area 
were  to  be  considered  from  zero  to  infinity,  might  be  2 to  2.5  times 
greater  (Reference  3).  However,  the  power  in  the  higher  frequencies 
contributes  very  little  to  the  response  of  an  aircraft. 

Power  spectra  measured  in  cumulus  clouds  are  shown  in  Figure  5, 
(Reference  2)  indicating  that  all  of  a wide  variety  of  rms  values 
demonstrate  the  characteristic  decrease  In  power  with  increasing  frequency. 
Analysis  of  numerous  spectra  under  varying  meteorological  conditions 
Indicates  that  a composite  of  all  possible  spectra  would  result  in 
a continuous  band  of  Intensities  varying  from  near  zero  to  possibly 
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as  high  as  16  feet  per  second  rms  (Reference  2).  Two  analytical 
representations  for  the  power  spectrum  of  atmospheric  turbulence  have 
received  extensive  use.  The  "Dryden"  spectrum 


I + 3L2fl.2 
(i+L2n2)2 


L = Scale  of  Turbulence 

Dryden  Spectra 


(14) 


i 


is  the  simplest  of  the  two  and  was  first  used  to  define  the  turbulence 
spectra  in  wind  tunnels.  The  von  Karman  spectra: 


<r*L  1 + ■§■ 

' [l  + (l.339La),']"/‘ 


05) 


von  Karman  Spectrum 


is  more  complicated,  but  is  a better  fit  to  the  recorded  spectra  for 
atmospheric  turbulence. 


\ 

( 

i 
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Figure  6 (Reference  2)  shows  recorded  spectra  and  spectra  that  have 
been  fitted  using  the  von  Karman  representation.  The  parameter  L is 
called  the  "scale  of  turbulence"  and  is  a measure  of  the  average  eddy 
size.  More  precisely,  t is  the  horizontal  separation  in  space  at  which 
the  gust  velocities  no  longer  display  any  correlation.  Considerable 
testing  of  the  von  Karman  spectrum  has  established  that  it  is  a 
reasonable  fit  to  all  levels  of  turbulence  and  that  the  value  of  L 
appears  to  be  a function  of  altitude.  Observations  indicate  that  the 
rms  turbulence  intensity,  a,  of  the  two  previously  mentioned  categories 
of  turbulence  also  vary  with  altitude  independent  of  each  other.  The 
first  category  is  referred  to  as  "non-storm"  and  the  second  category  is 
referred  to  as  "storm."  In  Figure  7 (Reference  4)  b-j  is  the  rms  turbulence 
intensity  averaged  over  the  time  that  the  aircraft  is  in  non-storm 
turbulence,  and  b2  is  the  corresponding  intensity  for  storm  turbulence. 
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As  indicated  by  b2,  storm  turbulence  is  most  severe  at  about  25,000  feet 
altitude  and  reaches  an  rms  intensity  level  of  almost  12  feet  per  second 
at  that  altitude.  The  probability  of  turbulence  (fraction  of  flight 
time  in  turbulence),  however,  is  greater  at  the  lower  altitudes.  In 
Figure  8 (Reference  4)  P-|  is  the  probability  of  non-storm  turbulence  at 
various  altitudes  and  P2  is  the  same  result  for  storm  turbulence. 
Especially  in  the  higher  altitudes,  the  greatest  portion  of  a flight 
will  be  in  smooth  air  as  indicated  by  the  fact  that  P1  + P2  < 1 . 


i 
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SECTION  III 
DYNAMIC  ANALYSIS 

When  a stable  system  Is  driven  by  a steady  sinusoidal  input,  the 
steady  state  response. will  be  a sinusoid  of  the  same  frequency,  and  the 
phase  angle  of  the  response  will  differ  from  that  of  the  excitation  by 
some  constant.  If  the  complex  amplitudes  of  the  excitation  and  response 
are  written  in  terms  of  the  steady  frequency  of  excitation,  then  the 
ratio  of  the  response  amplitude  to  the  excitation  amplitude  is  the 
"frequency  response  function."*  The  objective  of  the  dynamic  analysis 
is  to  formulate  the  frequency  response  function  for  a number  of  design 
parameters,  such  as  loads,  stresses,  or  accelerations  at  critical  locations 
on  the  aircraft. 

1 . OVERVIEW 

The  starting  point  for  a dynamic  analysis  is  a complete  description 
of  a model  of  the  aircraft  in  terms  of  geometry,  stiffness,  and  mass 
properties.  Normalized  free-free  airplane  modes  of  vibration  are  calculated 
together  with  their  associated  frequencies.  These  normal  modes  together 
with  rigid  body  modes  are  the  basis  for  the  generalized  coordinate 
system  of  the  analysis.  All  deflections  which  are  possible  for  this 
model  are  representable  as  a sum  of  the  products  of  the  mode  shapes,  and 
suitable  values  of  the  corresponding  generalized  coordinates.  The 
aerodynamic  surfaces  are  modeled  as  a series  of  "igid  panels,  each  of 
which  is  attached  to  the  structural  model.  As  the  structure  undergoes 
displacements,  the  panels  move  with  respect  to  the  free  airstream  and 
this  gives  rise  to  air  forces  on  the  various  panels.  The  generalized 


*The  "frequency  response  function,"  "transfer  function,"  and  "mechanical 
admittance"  are  terms  that  are  often  used  interchangeably. Strictly 
speaking,  the  frequency  response  function  and  the  mechanical  admittance 
refer  to  the  same  expression,  but  they  are  not  necessarily  the  same  as 
the  transfer  function.  The  transfer  function  will  be  the  same  expression 
only  for  the  special  case  of  a stable  linear/invariant  system  subjected 
to  steady  sinusoidal  excitation  and  viewed  after  all  starting  transients 
have  vanished.  However,  this  is  the  type  of  system  under  consideration. 
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mass  and  the  inertia  forces,  due  to  displacements  in  the  generalized  * 

coordinate  system,  are  calculated.  Because  they  are  frequency  dependent,  j 

generalized  response  and  excitation  air  forces  must  be  calculated  for  a ] 

distribution  of  frequencies  over  the  frequency  range  in  which  significant 
aircraft  response  is  possible.  j 

The  equations  of  motion  are  then  formulated  and  solved  for  each  j 

frequency  to  give  the  frequency  response  in  each  of  the  generalized  \ 

coordinates.  The  responses  are  the  independent  variables  in  a series  ] 

of  load  equations  that  yield  the  design  parameters  which  are  the  ob- 
jective of  the  dynamic  analysis.  ] 

i 

2.  LUMPED  PARAMETER  SYSTEM  j 

The  aircraft  structure  is  a continuous  distribution  of  mass  and  j 

stiffness,  and  calculations  of  normal  modes  and  frequencies  on  the  1 

basis  of  such  a system  is  a formidable  task.  The  equations  of  motion 
of  a continuous  system  involve  partial  differential  or  integral  equations  | 

which  are  relatively  difficult  to  solve.  Creation  of  a lumped  parameter  I 

mathematical  respresentation  of  the  structure,  or  the  use  of  finite  ] 

elements,  reduces  the  problem  to  one  involving  a finite  number  of  1 

simultaneous  ordinary  differential  equations,  which  are  easier  to  solve. 

The  basis  for  a lumped  parameter  or  finite  element  representation 
of  a continuous  system  is  the  assumption  that  deformations  of  the  continuous 
structure  may  be  approximated  by  deflections  and  rotations  at  a finite 
number  of  discrete  points,  and  that  the  inertial  properties  of  the 
structure  may  be  represented  by  suitable  masses  and  inertias  locatea  at 
these  points.  In  the  case  of  a wing  of  low  aspect  ratio,  a two-dimensional 
array  of  points  is  often  used  in  order  to  adequately  represent  the 
system.  A high  aspect  ratio  wing  which  is  very  stiff  in  the  chordwise 
direction  may  often  be  adequately  represented  by  points  on  a single 
line  located  along  the  elastic  axis,  with  concentrated  mass  and  inertia 
at  each  point. 
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3.  MASS  AND  STIFFNESS  MATRICES 

The  equations  of  motion  for  free  vibration  are  most  readily 
formulated  in  matrix  form  with  mass  and  stiffness  matrices.  The  mass 
matrix  for  a lumped  parameter  model  is  diagonal,  with  the  lumped  masses 
and  inertias  forming  the  diagonal.  The  masses  are  assumed  to  undergo 
deflections  in  only  one  direction  — normal  to  the  plane  of  the  surface, 
and  inertias  undergo  rotations  about  only  one  axis.  The  stiffness 
matrix  consists  of  both  diagonal  and  off-diagonal  terms  and  is  somewhat 
difficult  to  calculate  directly.  The  flexibility  matrix,  which  is  the 
inverse  of  the  stiffness  matrix,  is  relatively  easy  to  formulate 
analytically,  or  to  determine  experimentally  if  the  structure  is  available 
for  test.  The  i^1  column  of  the  flexibility  matrix  is  the  array  of 
deflections  (or  rotations)  at  each  of  the  mass  points  (or  inertia  locations), 
due  to  a unit  load  (or  moment)  applied  to  the  i^1  mass  point  (or  Inertia 
location).  Having  formulated  the  flexibility  matrix,  the  stiffness 
matrix  may  be  obtained  by  matrix  inversion.  The  finite  element  method 
allows  the  stiffness  matrix  to  be  formulated  directly,  and  is  often  used 
for  large  structures. 

4.  NORMAL  MODES 

If  structural  damping,  which  has  relatively  small  effects  on 
aircraft  mode  shapes  and  frequencies  is  neglected,  the  differential 
equation  for  free  vibration  of  an  aircraft  structure  may  be  written  as: 


[“<]{*}  + [“*]{«}  '{»}  06 

where 

= Mass  matrix 
= Stiffness  matrix 

The  assumption  that  the  response  of  the  system  is  simple  harmonic 
motion  reduces  this  differential  equation  to  the  following  algebraic 
equation  which  has  a non-trivial  solution  if,  and  only  if,  the  deter- 
minant of  the  coefficient  is  zero: 


[W-'NjW  * {°}  (17) 
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The  solution  to  this  equation  is  a set  of  natural  frequencies  , 
and  corresponding  eigenvectors  which  satisfy  the  above  equation.1 

When  a large  number  of  masses  are  used  to  represent  the  structure, 
the  problem  of  calculating  free-free  normal  modes  directly  becomes 
difficult.  In  this  case,  it  is  often  advantageous  to  calculate  modes 
and  frequencies  for  each  of  the  major  components  of  the  structure  individually, 
and  then  to  combine  these  so-called  "cantilever  modes"  together  with 
rigid  body  modes  in  a second  eigenvalue  problem  to  obtain  the  free- 
free  aircraft  modes.  In  deriving  the  cantilever  modes,  the  point  of 
attachment  of  the  particular  component  to  the  main  structure  is  assumed 
to  be  rigid.  It  is  important  that  all  of  the  cantilever  modes  whose 
frequencies  are  within  the  frequency  range  in  which  significant  aircraft 
response  is  possible  be  determined  in  order  to  preserve  accuracy  in  the 
final  free- free  modes. 


5.  EQUATIONS  OF  MOTION 

The  equations  of  motion  written  in  terms  of  the  generalized  co- 
ordinates for  the  aircraft  are  a set  of  linear,  nonhomogeneous,  sim- 
ultaneous differential  equations  of  sec""*  order.  The  coefficient 
matrices  consist  of  both  aerodynamic  and  structural  terms,  and  may 
be  written  as  follows: 

[ms+  *]{,}♦[*♦  ■>*]{«}  +[KS  + '<a]{i}-K}V»  08) 


where 


M-M-M 


= acceleration,  velocity,  and  displacement 
matrices  of  the  generalized  coordinates. 


gust  velocity 

modal  force  coefficient  matrix  per  unit  gust  velocity 
generalized  structural  mass  matrix 
generalized  aerodynamic  mass  matrix 
generalized  structural  damping  matrix 
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= generalized  aerodynamic  damping  matrix 
= generalized  structural  stiffness  matrix 
= generalized  aerodynamic  stiffness  matrix 

The  structural  coefficient  matrices  are  all  diagonal  matrices  as  a 
result  of  selecting  normal  modes  as  a basis  for  the  generalized  co- 
ordinate system.  The  generalized  mass  may  be  calculated  from  the 
following  summation  over  N masses  representing  the  entire  aircraft: 


iTi^i 


(19) 


« 
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where 

jL j = the  deflection  at  the  lumped  mass  due  to  a unit 
deflection  in  the  1*h  generalized  coordinate 

- mass  of  the  jtl)  lumped  mass 

For  the  rigid  body  translation  modes,  0^  = 1,  and  the  resulting  term 
in  the  generalized  mass  is  equal  to' the  mass  of  the  aircraft  The 
coordinate  system  for  rigid  body  rotation  modes  must  be  the  principal 
axis  of  the  aircraft  in  order  to  preserve  orthogonality.  For  these 
modes,  the  generalized  mass  term  is  the  normalized  mass  moment  of 
inertia  about  the  appropriate  principal  axis. 


The  structural  contribution  to  the  generalized  stiffness  Is  cal- 
culated from  the  generalized  mass  as  follows: 


Ks  = Mg  <*>* 
*11  Nl 


(20) 


where 

Uu  = natural  frequency  of  the  itfl  normal  mode. 

N1 

The  natural  frequency  for  the  rigid  body  modes  is  zero,  hence  the 
generalized  stiffness  matrix  is  singular. 
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Stri'ctural  damping  results  from  friction  within  the  structure, 
particularly  at  joints  between  components.  Unlike  viscous  damping 
forces,  these  damping  forces  are  proportional  to  the  strain  of  the 
system,  yet  are  in  phase  with  the  velocity.  This  contribution  may 
be  formulated  as  follows: 

Ml*}* ' [•][«•]{«}. 

As  a result  of  assuming  simple  harmonic  motion: 

{’}  *'"{«} 

Substitution  gives: 

Ml*}  ■ -A- Mis] 

Hence: 

Is]  ■ -A- Mis] 

The  damping  coefficient  matrix  £gj  is  determined  from  vibration 
tests,  when  available;  or  is  assumed  to  be  about  5%  of  critical  damping. 

6.  LOAD  EQUATIONS 

Load  equations  for  calculating  the  design  parameters  which  are 
the  objective  of  the  dynamic  analysis  are  similar  in  form  to  the 
equations  of  motion,  and  may  be  expressed  as: 


(21) 


(22) 


(23) 


(24) 


{«}  s [mr]{<>}  ♦[*]{*}♦[*]{«  } ♦ {Cr}  "g  M (25) 

where 

Mr  = generalized  mass  matrix  for  response  load  equations 
Dr  = generalized  damping  matrix  for  response  load  equations 
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[l<r]  = generalized  stiffness  matrix  for  response  load  equations 

{c  l = modal  force  coefficient  matrix  for  response  load 
*■  rJ  equations 

These  equations  are  formulated  by  transforming  the  generalized  inertia, 
damping,  and  stiffness  forces  Into  loads,  stresses,  accelerations, 
etc.,  per  unit  deflection  of  the  generalized  coordinates. 

7.  AERODYNAMIC  FORCES 

The  generalized  external  forces  which  result  from  the  motion  of 
the  system  may  be  expressed  as  a column  matrix  of  modal  forces  {qd{ 
These  forces 

1 
j 

i 

where 

‘ M- 

M- 

H- 

Selecting  an 

generalized  external  force  results  in  a matrix  with  complex  coeffi-’ 
dents.  The  assumption  of  simple  harmonic  motion  allows  this  expression 
to  be  successively  expanded  as: 


■ [•]{<}  + •£■[►]{*}  (27) 


may  be  calculated  from  the  following  matrix  equation: 

(26) 

transformation  from  displacements  of  the  generalized 
coordinate  system  to  displacements  at  the  aerodynamic 
panels,  generally  referred  to  as  mode  shapes. 

complex  transformation  from  two-dimensional  panel  forces 
to  three-dimensional  panel  forces  which  includes  finite 
time  delays  for  the  transmittal  of  induced  lifts, 
generally  referred  to  as  dynamic  induction  matrix. 

oscillatory  panel  forces  per  unit  sinusoidal  panel 
displacements,  generally  referred  to  as  lift  and  moment 
matrix. 

excitation  frequency,  u,  and  evaluating  the  foregoing 
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The  aerodynamic  stiffness  matrix, 
which  are  in  phase  with  the  genera 
dynamic  damping  matrix,  , results  from  those  air  forces  which  are 
90-degrees  out  of  phase  with  the  generalized  coordinates.  Hence: 

[ka  ] = REAL  PART  OF  (28) 

and: 

[da]  = IMAGINARY  PART  OF  [$T][g  ][  A ][  + ] (29) 


«aJ  , results  from  those  air  forces 
ized  coordinates,  and  the  aero- 


A1 though  the  mass  of  the  air  in  close  proximity  to  the  aerodynamic 
surface  does  impart  inertial  forces  to  the  surface  itself,  this  con- 
tribution is  negligible  compared  with  the  mass  of  the  structure. 
Therefore,  is  assumed  to  be  zero: 

[«.]■[»]  <»> 


The  remaining  external  forces  arise  from  the  gust  induced  angle  of 
attack  ag(t).  A phase  angle  delay,  which  is  proportional  to  the 
panel  location  aft  of  some  reference  point,  appears  in  the  expression 
in  order  that  the  effects  of  the  gust  will  be  felt  progressively 
along  the  flight  path  at  a speed  equal  to  the  velocity  of  the  air- 
craft. Hence,  for  a sinusoidal  gust  field  of  frequency  w and  maximum 
amplitude  wg,  the  gust  induced  angle  of  attack  for  the  various  panels 
may  be  expressed  as: 


iwt 


{•*} 


(31) 


where 

Vj  = true  velocity  of  the  aircraft 

9 = rr*  - phase  angle  delay  for  each  panel 
VT 

A = panel  location  aft  of  some  aircraft  reference  point 
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The  form  of  the  final  equation  for  gust  induced  external  forces, 
is  similar  in  form  to  the  expression  for  the  response  external  forces, 
the  essential  difference  being  the  [d]  {q}  is  replaced  by  {ag(t)}. 
Hence: 

{«.}  * [*T][S  ][*•]{<■«">}  <32) 


where 

[Ag]  = oscillatory  panel  forces  per  unit  sinusoidal  angle 
of  attack 

From  this  equation,  the  complex  column  matrix  {C}  in  the  equations 
of  motion  is  given  by: 

{ = }»  -J-  [*T][s][*, ]{.'*}  (33) 

8.  AERODYNAMIC  THEORIES 

A number  of  aerodynamic  theories  of  various  degrees  of  complexity 
are  available  to  calculate  the  oscillatory  panel  forces.  The  "doublet- 
lattice"  (Reference  5)  is  an  example  of  one  of  the  more  complex  theories, 
and  calculations  based  on  this  theory  agree  well  with  tests.  The 
aerodynamic  surfaces  are  divided  into  a series  of  trapezoidal  "panels" 
by  a series  of  lines  parallel  to  the  free  stream,  and  a series  of 
lines  of  roughly  constant  percent  chord.  On  each  panel,  a steady 
"horseshoe  vortex"  is  used  to  develop  the  steady  portion  of  the  lift, 
and  a line  of  oscillatory  "doublets"  is  used  to  represent  the  oscillatory 
portion.  The  horseshoe  vortex  is  located  so  that  the  bound  segment  of 
the  vortex  coincides  with  the  quarter-chord  of  the  panel  and  the  parallel 
segments  of  the  vortex  trails  from  the  ends  of  the  quarter-chord  of  the 
panel  in  the  free  stream  direction.  The  lines  of  oscillatory  doublets 
coincide  with  the  quarter-chord  of  the  panel.  "Downwash  collocation" 
points  are  located  on  the  three-quarter-chord  of  each  panel,  centered 
spanwise.  The  downwash  from  each  of  the  vortices  and  doublets  is  summed 
at  each  collocation  point  and  the  component  normal  to  the  surface  is 
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equated  to  zero.  This  results  in  a system  of  simultaneous  equations,  the 
solution  of  which  gives  the  aerodynamic  forces  on  each  of  the  panels, 
per  unit  oscillatory  deflection  at  a given  frequency. 

9.  SOLUTION  TO  THE  EQUATIONS  OF  MOTION 

A solution  for  the  equations  of  motion  Is  required  for  a distribution 
of  frequencies  from  near  zero  to  the  highest  frequency  at  which  significant 
structural  response  is  possible.  Combining  the  structural  and  aerodynamic 
terms  into  single  unsubscripted  terms  yields  the  following  equation: 

["]{«} + [D]{*}  + [K]{<,}s‘{co}*«(,)  <34) 

The  response  of  a linear  system  to  steady  sinusoidal  excitation  will 
be  simple  harmonic  motion  at  the  same  frequency  as  the  excitation, 
but  with  some  finite  phase  angle  between  the  excitation  and  response. 

Hence: 

{«}-k}*“”+*  (35) 

when 

WgiDse1"1  (36) 

where  {qQ>  is  the  maximum  amplitude  of  the  steady  sinusoidal  response, 
and  <p  is  the  phase  angle  between  the  excitation  and  response.  Successive 
differentiation  of  the  modal  displacement  gives  the  modal  velocity  and 
acceleration  as: 

+*  (17) 

and 

{q}  * - «2{<0}  +*  (38) 
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{ 

i 

\ 

i 


Substituting  the  modal  displacement,  velocity,  and  acceleration  into  the 
equation  of  motion  and  dividing  by  elwt  + ^ yields: 

[tK]  " w2[M]  + ,w[D]]{s}*-{co}  (39) 


From  which  the  solution  is: 


<%}  ■-[[«]  --•[-]  * -wJk} 

The  right-hand  side  of  the  foregoing  equation  is  the  frequency  response 
function  of  that  part  of  the  system  which  relates  the  modal  displacement 
to  the  sinusoidal  gust  velocity: 

{hOw)}  = £[k]  +0«*i)2[m  ] - (itt»)[l>]J{c0]'  (41) 


Frequency  Response  Function 


If  it  were  not  for  the  aerodynamic  contribution  to  the  damping  and 
stiffness  matrices,  the  frequency  response  function  for  this  system 
would  display  a single  peak  as  does  the  frequency  response  function 
for  a single  degree-of-freedom  damped  oscillator,  which  when  normal 
ized  appears  as  in  Figure  9 (Reference  6). 


Figure  9.  Frequency  Response  Function  of  Simple 
Damped  Oscl’lator 


28 


AFFDL-TR-76-162 


where 

H (iu>)  = frequency  response  function  for  single  degree-of-freedom 
damped  oscillator 

? = ratio  of  actual  damping  to  critical  damping 
<*)n  = undamped  natural  frequency 
co  = sinusoidal  excitation  frequency 

The  significance  of  the  off-diagonal  contribution  of  the  aerodynamic 
terms  to  the  damping  and  stiffness  matrices  is  to  create  cross- 
coupling between  the  normal  modes.  Figure  10  (Reference  7)  shows  the 
frequency  response  function  of  the  modal  displacement  of  the  first  sym- 
metrical flexible  free-free  mode  of  vibration  of  a conmercial  transport 
aircraft. 

i 

i 

i 10.  SOLUTION  TO  THE  LOAD  EQUATIONS 

A solution  to  the  load  equations  Is  obtained  from  the  solution 
' to  the  equations  of  motion,  using  the  same  assumptions  of  simple  harmonic 

motion.  Substituting  the  modal  displacement,  velocity,  and  acceleration 
written  in  terms  of  {qQ>  into  the  load  equations,  collecting  terms, 

and  dividing  by  eiu)t  + ^ yields: 

k}  *[[kr] + ]]{«•«}  +{cr}  <«> 

where 

{Rq}  is  the  maximum  amplitude  of  the  steady  sinusoidal  response 
and 

{qQ}  is  the  solution  to  the  equations  of  motion. 

The  frequency  response  functions  for  the  design  parameters  display 
the  same  multi-peak  characteristics  as  the  frequency  response  function 
of  the  normal  coordinates.  Figure  11  shows  the  frequency  response 
function  for  wing  bending  moment  at  1/3  span  of  a commercial  transport 
aircraft  subject  to  vertical  gust  excitation  of  one  foot  per  second 
rms. 
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Figure  10.  Frequency  Response  Function  of  First  Free-free 
Airplane  Mode  of  Vibration 


ft  (rad/ft) 

Figure  11.  Frequency  Response  of  Wing  Bending 
Moment  at  1/3  Span 
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SECTION  IV 
SPECTRAL  ANALYSIS 

For  a linear  system  excited  by  a random  function  which  is  both 
Gaussian  and  stationary,  the  PSD  of  the  output  is  expressed  in  terms 
of  the  PSD  of  the  input  by  the  following  equation: 

<t>  = (43> 

output  1 1 Input 

where 

H = frequency  response  function 

The  relative  simplicity  of  this  equation  has  motivated  the  assumptions 
that  atmospheric  turbulence  may  be  characterized  as  a stationary 
random  process  whose  gust  amplitudes  have  a Gaussian  distribution,  and 
that  the  flight  of  an  aircraft  through  the  atmosphere  may  be  characterized 
as  a linear  system  represented  by  a frequency  response  function. 

The  concept  of  stationarity  that  is  implied  is  that  the  statistical 
properties  of  the  input  averaged  over  a short  time  interval  do  not  vary 
significantly  from  one  interval  to  another.  This  concept  is  sometimes 
called  "self-stationarity"  to  avoid  confusion  with  the  more  classical 
definition  of  stationarity  which  involves  averages  over  an  "ensemble" 
(Reference  8)  of  functions,  rather  than  averages  in  time  along  a single 
function. 


If  X(t)  characterizes  a velocity  component  of  atmopsheric  turbulence, 
then  the  Gaussian  assumption  implies  that  the  probability  distribution 
of  X(t)  may  be  expressed  in  terms  of  two  parameters  --  the  "variance" 
and  the  "mean  value."  This  distribution  has  a bell  shape  given  by  the 
following  expression: 
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where 


T 


f X(t)  dt  = mean  value  of  X(t) 

T o 

T 2 

dt  = variance  of  X(t) 


o 


p(x)  = probability  density  of  X(t) 


Accepting  these  assumptions,  the  power  spectral  relation  for  the 
response  of  an  aircraft  becomes: 

i ,2 

i («*»)  = |h  (iw)| 

R 

where 

<P (co)  = PSD  of  the  atmosphere  as  detailed  in  Section  II 

H(iu))  = frequency  response  function  of  the  design  parameters 
as  detailed  in  Section  III 

4>r(oj)  = PSD  of  the  design  parameter 


Two  real  and  positive  parameters  are  calculated  from  the  Dower  spectrum 
of  each  of  the  design  parameters  for  use  in  the  design  process.  The 
first  parameter.  A,  is  the  rms  response  per  unit  rms  turbulence  intensity, 
and  it  is  the  square  root  of  the  area  under  the  response  parameter 
power  spectrum.  The  second  parameter,  Nq,  is  a characteristic  frequency, 
calculated  from  the  first  moment  of  the  area  under  the  PSD  curve  about 
the  zero  frequency  line.  It  is  the  expected  number  of  exceedances  of 
the  zero  load  level  per  unit  time.  These  expressions  are: 


and: 


A 


d>  (o>)  dui 
R 


"o' 


4 


cu  d)  ((u)dtoi 
R 


2 ir  A 
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1.  DESIGN  ENVELOPE  ANALYSIS 

In  the  PSD  design  envelope  analysis,  it  is  assumed  that  the  ratio 
of  the  "correct"  design  value  to  the  rms  design  value  is  a constant  n^. 

If  the  design  value  for  rms  turbulence  intensity  is  ow,  then  the  rms 
response  is  aJS,  so  that  the  "correct"  design  value  is  (Reference  9): 

(48) 

In  practice,  ow  and  need  not  be  considered  separately  and  their 
product  is  defined  as  Ua,  the  nominal  turbulence  intensity.  The  design 
value  for  the  design  parameters  becomes: 

yd  « Ua A (49) 

Acceptable  values  of  may  be  determined  from  the  analysis  of  previous 
aircraft  which  were  gust  critical,  and  have  demonstrated  a long  service 
life.  A is  calculated  and  the  foregoing  equation  is  solved  for  Uo  when 
^design  is  set  equa1  t0  the  allowable  value  of  the  design  parameter: 

y , yi*~  (50) 

* A 

where 

yLA  = allowable  value  for  design  parameter 

y-|_g  = value  of  the  design  parameter  in  steady  level  (1-g)  flight 

A - rms  response  of  design  parameter  per  unit  rms  gust 
excitation 

In  an  effort  to  establish  acceptable  values  of  Ua,  the  Federal  Aviation 
Adnlnlstratlon  sponsored  a study  of  a number  of  proven  aircraft  designs. 
Company  engineering  staffs  determined  the  value  of  Ua  that  would 
produce  the  limit  allowable  stress  value  at  the  most  critical  location 
on  each  aircraft.  This  amounts  to  selecting  the  lowest  value  of  U0 
from  the  analysis  of  all  of  the  locations  which  were  considered  to  be 


i : 

J 
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potentially  critical.  The  results  of  two  such  analyses  are  listed 
below  (Reference  4): 


AIRCRAFT  MODEL  Uo  VERTICAL  GUST  ALTITUDE 


Boeing  707-720B  111  ft/sec 

Lockheed  749  88  ft/sec 


22000  ft 
7000  ft 


i 

j 

■1 

i 


Based  on  this  type  of  analysis,  the  following  values  of  Ua  as  a function 
of  altitude  has  been  established.  The  solid  lines  are  proposed  standards 
applicable  to  commercial  transport  category  aircraft,  and  the  dotted 
line  is  the  requirement  of  the  US  Air  Force  (Reference  10). 


Figure  12.  Derived  Gust  Velocity  for  Design 
Envelope  Analysis 


where 

Vfi  = gust  penetration  velocity  (civil  aircraft) 

Vu  = maximum  continuous  velocity  in  steady  level  flight 
" (military  aircraft) 

Vc  = maximum  cruise  velocity  (civil  aircraft) 

Vq  = maximum  dive  velocity  (civil  aircraft) 
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2.  MISSION  ANALYSIS 

The  mission  analysis  permits  the  intended  usage  of  the  aircraft  to 
be  reflected  in  the  design.  The  objective  of  the  mission  analysis  is 
to  establish  design  levels  so  that  during  the  intended  operation  of  the 
aircraft,  the  expected  number  of  exceedances  of  specific  load  levels 
per  hour  of  operation  do  not  exceed  an  acceptable  value. 

The  basis  for  the  mission  analysis  is  the  assumption  that  both 
storm  and  non-storm  turbulence  occur  in  patches,  where  each  patch  is 
a stationary,  Gaussian  random  process  with  a particular  rms  turbulence 
intensity  o.  Although  a is  a constant  within  each  patch,  it  is  assumed 
that  the  distribution  of  o's  from  the  various  patches  in  both  categories 
of  turbulence  is  Gaussian.  For  a single  patch,  the  number  of  exceedances 
per  unit  time  of  level  y Is  given  by  Equation  51.  This  equation  is 
derived  in  Appendix  A using  Houbolt's  approach  (Reference  11). 


where 


N(y) 


N = expected  number  of  exceedances  of  the  level  y = 0 per 
unit  of  time 


(51) 


0^  = rms  response  of  design  parameter 

y = level  of  response  design  parameter  for  which  exceedances 
are  calculated 


Invoking  the  assumption  that  the  various  values  of  a from  each  category 
have  a Guassian  distribution  results  in  the  following  modification  of 
the  foregoing  equation: 


where 


A 


N(y)  = N0 


♦ P2 


(52) 


probability  of  encountering  turbulence  of  non-storm  and 
storm  category  respectively 

rms  turbulence  intensity  of  non-storm  and  storm  turbulence 
respectively 

rms  response  of  the  design  parameter  per  unit  rms 
turbulence  excitation 
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A detailed  derivation  of  the  foregoing  equation  is  presented  in  Appendix  A. 
In  the  mission  analysis  design  process.  It  is  assumed  that  the  operation 
of  the  aircraft  consists  of  a specified  number  of  flights  where  each 
flight  is  one  of  a limited  number  of  "missions."  The  ratio  of  the  number 
of  flights  of  the  various  missions  to  the  total  number  of  flights  is 
called  the  "mission  mix,"  which  is  specified  as  part  of  the  design 
requirements.  Examples  of  some  missions  that  might  be  included  in  the 
design  of  fighter  aircraft  include: 

1 . Ferry 

2.  Training 

3.  Combat 

4.  Reconnaissance 

Each  mission  is  created  by  linking  together  a number  of  “mission  seg- 
ments" which  is  often  illustrated  in  a plot  of  altitude  vs  time  such 
as  the  plot  on  the  following  page.  Typical  mission  segments  which  might 
be  included  in  one  or  more  missions  are: 

1 . Takeoff 

2.  Climb 

3.  Cruise 

4.  Loiter 

5.  Air  to  air  combat  maneuvers 

6.  Air  to  ground  combat  maneuvers 

7.  Descent 

8.  Landing 

Having  decided  upon  how  many  flights  of  each  mission  will  be  flown  during 
the  life  of  the  aircraft,  and  how  each  mission  is  assembled  in  terms 
of  mission  segments,  a composite  mission  is  established  which  has  the 
same  ratio  of  the  various  mission  segments  as  the  sum  of  all  of  the 
flights  planned  for  the  lifetime  of  the  aircraft.  If  T is  the  duration 
of  the  composite  mission  and  T..  is  the  duration  of  the  ith  mission 
segment  in  the  composite  mission,  then  define  t.  to  be  that  fraction  of 

t h 1 

time  for  the  i n mission  segment  in  the  composite  mission: 

T, 

»,  8 -f-  (53) 
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Figure  13.  Mission  Profile  of  Training  Mission 


The  average  number  of  exceedances  per  unit  time  for  the  composite 
mission  will  be  given  by: 

Mrt-  IN.J,,  P,.^]  (54> 

mission  1 
segments 

A plot  of  N(y)  vs  y results  In  an  exceedance  curve  of  the  following  form: 


linaar  scale 
7 

Figure  14.  Typical  Exceedance  Curve 
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There  is  no  practical  load  level  for  which  the  number  of  exceedances 
will  be  zero,  hence,  in  order  to  utilize  the  above  result,  an  accep- 
table number  of  exceedances  of  some  particular  load  level  must  be 
established. 

Company  engineering  staffs  have  calculated  the  expected  number  of 
exceedances  of  limit  strength  per  hour  of  operation  of  a number  of 
successful  aircraft,  and  the  results  of  two  such  analyses  are  as 
follows  (Reference  4): 

EXCEEDANCES  OF  LIMIT  LOAD  PER  HOUR 
AIRPLANE  MODEL  Vertical  Gust  Lateral  Gust 

Boeing  707-720B  1.1  x 10~5  4.0  x 10‘6 

Lockheed  749  1.8  x 10~5  2.5  x 10“4 

The  appropriate  design  level  for  the  nuuber  of  exceedances  of  limit 
load  per  hour  of  operation  due  to  vertical  gust  was  selected  to  be 
(Reference  10): 

N(y)  = 2.0  « IO’S  (55) 

which  amounts  to  one  exceedance  of  limit  load  per  50,000  hours  of 
operation.  This  is  currently  the  requirement  for  US  Air  Force  aircraft 
and  is  the  proposed  standard  for  civil  transport  aircraft. 
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SECTION  V 

NONSTATIONARY  TURBULENCE  MODELS 

The  stationary  assumption  has  been  made  in  order  to  simplify 
aircraft  response  calculations,  but  this  assumption  is  not  entirely 
satisfactory.  The  initial  encounter  with  a gust  field  is  often  char- 
acterized by  rapid  changes  in  the  statistical  properties  of  the  aircraft 
response  parameters  with  time,  and  these  properties  may  vary  continuously 
while  the  aircraft  is  in  turbulence.  This  effect  cannot  be  accounted 
for  in  a stationary  analysis.  Indeed,  the  analysis  of  real  data  indicates 
that  the  traditional  PSD  procedure  consistently  underestimates  the  number 
of  exceedances  of  the  higher  load  levels  due  to  turbulence., 

The  degree  to  which  the  traditional  analysis  underestimates  the 
number  of  exceedances  in  these  higher  load  levels  has,  on  occasion, 
been  exaggerated.  Some  of  the  reports  on  the  inadequacy  of  the  traditional 
analysis  are  the  result  of  a lack  of  understanding,  and  it  is  appropriate 
to  expose  some  of  these  misconceptions. 

The  most  flagrant  misconception  is  that  the  present  analysis  is 
entirely  Gaussian.  The  number  of  exceedances  versus  load  levels  for 
a single,  stationary,  Gaussian  process  is  compared  to  points  indicative 
of  real  data  in  Figure  15. 


Figure  15.  Gaussian  Exceedance  Curve 
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The  single  stationary,  Gaussian  process  plots  as  a straight  line 
in  Figure  15,  whereas  a distribution  of  Gaussian,  stationary  patches 
plots  as  shown  by  the  dashed  line,  resulting  in  a considerably  better 
fit  to  the  points  indicative  of  real  data. 

On  other  occasions,  the  fact  that  the  Gaussian  Patch  model  assumes 
that  the  mean  square  intensity  of  the  various  patches  come  from  two 
different  processes  (storm  and  non-storm),  has  not  been  recognized. 

The  significance  of  this  assumption  appears  in  Figure  16. 


Figure  16.  Gaussian-Patch  Exceedance  Curve 

In  this  figure,  the  Gaussian-patch  model  from  a single  distribution 
is  the  solid  straight  line.  However,  when  two  distributions  are 
considered,  it  plots  as  the  dashed  line.  Although  the  use  of  two  dis- 
tributions give  results  which  are  considerably  better,  the  number  of 
exceedances  of  the  higher  load  levels  is  still  underestimated.  However, 
the  difference  is  usually  considered  to  be  negligible. 

Other  analyses  of  data  have  permitted  three  distributions  to  be 
used,  where  the  third  distribution  accounts  for  the  very  rare  occurrence 
of  extreme  turbulence  (Reference  12).  In  view  of  the  added  complexity 
of  including  a third  distribution,  the  value  of  such  a distribution  can 
be  debated. 
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1.  TRANSIENT  OVERLOADS 

In  the  Gaussian-patch  analysis,  the  transients  resulting  from  the 
transition  from  one  patch  to  another  are  not  considered.  The  stationarity 
assumption  has  been  questioned  on  the  basis  that  it  may  ignore 
significant  transient  overloads.  At  least  one  study  has  indicated  that 
transient  response  in  excess  of  30%  above  the  stationary  response  has 
been  observed  (Reference  13).  The  analysis  utilized  a nonstationary 
model  for  turbulence  of  the  following  form: 

w(t)  * o(t)  z ( t ) (56) 

where 

a(t)  = a deterministic  modulating  function 

z(t)  = spectrally  shaped,  Gaussian  white  noise 

w(t)  = nonstationary,  possibly  non-Guassian  turbulence 

The  modulating  function  a(t)  that  was  used  was  a ramp  function  of  the 
following  form: 


Figure  17.  Modulating  Function 


The  high  transient  was  observed  utilizing  an  enveloping  function  with 
an  extremely  short  rise  time  (90%  in  .01  seconds)  as  compared  with  the 
natural  response  time  of  the  aircraft  model.  In  effect,  the  aircraft 
was  probably  subjected  to  a severe  discontinuity  as  indicated  in 
Figure  18. 
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Figure  18.  Nonstationary  Excitation  Function 


The  dashed  line  represents  the  function  z(t)  and  the  solid  line  represents 
the  product  of  z(t)  and  a(t). 

In  summary,  the  Gaussian-patch  model  fails  to  predict  the  number 
of  exceedances  of  the  higher  load  levels,  but  the  magnitude  of  this 
deficiency  is  small  for  many  applications  to  gust  design  problems. 

2.  UNIFORMLY  MODULATED  PROCESS 

Past  efforts  to  model  atmospheric  turbulence  as  a nonstationary 
process  have  included  the  notion  that  turbulence  may  be  represented 
as  a stationary  Gaussian  random  process  modulated  by  an  enveloping 
function  as  follows: 


w(t)  = f(t)  z(t) 


(57) 


where 

f ( t ) = discrete  modulating  function 

z(t)  = Gaussian  white  noise  that  has  been  spectrally  shaped 
w(t)  = nonstationary,  possibly  non-Gaussian  turbulence 

A variety  of  enveloping  functions  can  be  devised  which  are  thought 
to  represent  a number  of  different  types  of  nonstationary  turbulence 
encounters.  This  model  is  referred  to  as  "uniformly  modulated"  and 
is  particulary  attractive  because  a Fourier  series  may  be  used  to 
express  the  enveloping  function,  and  equations  for  the  statistical 
properties  of  the  response  parameters  (including  the  number  of  ex- 
ceedances per  unit  time)  may  be  derived  in  closed  form. 
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Another  approach,  which  could  be  considered  to  be  an  expansion 
of  the  uniformly  modulated  process,  is  to  model  turbulence  in  the 
following  form  (Reference  14): 

w(0  = f(l)  1(1)  + 9(1)  (58) 

where  z,  f,  and  g are  independent  stationary  Gaussian  processes.  In 
effect,  g(t)  is  the  time  varying  mean  value  of  the  turbulence  and 
f(t),  although  random,  performs  a similar  function  to  that  of  f(t) 
in  the  uniformly  modulated  process;  z(t)  is  the  same  for  both  models. 

This  more  complicated  model  has  an  advantage  in  that  it  gives  a less 
uniform  structure,  and  thus  has  an  "appearance"  more  like  that  of  real 
turbulence. 

3.  TIME  VARYING  POWER  SPECTRA 

Other  approaches  which  are  computationally  more  difficult,  embody 
the  concept  of  a time-varying  power  spectrum.  The  "evolutionary  power 
spectrum"  first  presented  by  M.  B.  Priestly  is  of  the  form  (Reference  15): 

* |A(1,X1)|2  «#»(m  <59> 

where 

d(fi)  = stationary  power  spectrum  of  the  atmosphere,  i.e.,  the 
von  Karman  or  Dryden  spectrum 

A(t,n)  = time  varying  transfer  function 

d(t,fi)  = evolutionary  power  spectrum 

The  time-varying  transfer  function  performs  a similar  operation  in 
the  frequency  domain  to  what  the  enveloping  function  of  the  uniformly 
modulated  function  performed  in  the  time  domain. 

W.  D.  Mark  has  made  use  of  the  fact  that  the  autocorrelation  function 
and  the  PSD  function  are  Fourier  transform  pairs  (Reference  16).  For 
the  stationary  case; 

R(t)  s J $fa)4iu 
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and: 

#“»>  ’ ~2w  /#',WT  dT  (61) 

o 

Mark  obtained  a time  varying  autocorrelation  function  and  then  utilized 
the  Fast  Fourier  Transform  (FFT)  technique  to  obtain  the  time  varying 
"instantaneous  power  spectrum." 

Recently,  Mark  and  Fisher  investigated  the  effects  of  nonstationary 
behavior  on  the  spectra  of  atmopsheric  turbulence  using  the  instantaneous 
power  spectrum  (Reference  16).  They  concluded  that  the  shape  of  the 
von  Karman  spectrum  will  be  virtually  unchanged  when  the  enveToping 
function  of  the  uniformly  modulated  process  is  sufficiently  well 
behaved  that: 


(l«  <r{ t) ) j < 0.04 


(62) 


where 

L = scale  of  turbulence  of  the  von  Karman  spectrum 
o(t)  = time-varying  mean  value  (enveloping  function) 

Mark  and  Fisher  arrived  at  this  conclusion  by  starting  with  an  ex- 
pression for  the  instantaneous  power  spectra  of  a uniformly  modulated 
process: 

t)  = [<r(t)j  $(&)  (63) 

They  derived  a new  series  expansion  for  the  instantaneous  spectra  and 
wrote  the  first  two  terms  as: 

♦ (a.i)«o^t)[4(0)--jj^i--^-  (u.<M0)**(O)]  (64) 

Note,  that  this  results  in  the  regular  power  spectra  when  o(t)  is  con- 
stant, and  thus  the  second  term  displays  the  results  of  any  non- 
stationary fluctuations  in  the  modulating  function  o(t). 


j 

i 

j 

1 


j 

I' 

4 


j 
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In  order  to  determine  the  magnitude  of  the  second  term  necessary 
to  produce  a noticeable  change  in  the  power  spectrum,  Mark  and  Fisher 
generated  a number  of  nonstationary  samples  (i.e.,  time  series)  from 
a uniformly  modulated  process,  and  compared  the  power  spectra  of  these 
samples  to  the  power  spectra  of  "Gaussian  equivalent"  samples.  In  this 
investigation  Gaussian  equivalent  samples  were  constructed  such  that  the 
crossings  of  the  zero  load  level  in  time  were  identical  to  that  of  the 
nonstationary  samples,  yet  the  distribution  of  amplitudes  was  Gaussian. 
Mark  and  Fisher  utilized  the  so-called  "arc-sine  law"  to  construct  the 
Gaussian-equivalent  samples.  This  law  states  that: 


R,(r)  * sin[-|-  r0j(t)] 


where 


R (t)  = autocorrelation  function  of  a Gaussian  (possibly  non- 
stationary) random  process  z(t) 


(x)  = autocorrelation  function  of  a process  whose  magnitude  is 
+1  when  z(t)  > 0,  and  -1  when  z(t)  < 0.  Therefore,  the 
zero  crossings  will  be  identical  to  z(t). 


Figure  19  illustrates  the  procedure  used.  The  magnitude  of  the  first 
term  was  compared  to  the  value  of  the  second  term  in  making  the 
evaluation. 

4.  SELF-SIMILARITY 

Self-similarity  is  a concept  which  may  some  day  produce  a model 
which  characterizes  both  the  nonstationary  and  non-Gaussian  aspects 
of  atmospheric  turbulence.  Whereas,  the  concept  of  stationary  implies 
that  the  turbulence  is  invariant  with  respect  to  addition  of  time, 
the  concept  of  self-similarity  Implies  that  turbulence  is  invariant  with 
respect  to  multiplication  of  time  (Reference  17).  Specifically,  a 
process  with  random  variable  x is  self-similar  if  there  exists  some 
function  F(h)  such  that  the  probabilistic  structure  of  x is  unchanged 
when  the  wavelength  is  multiplied  by  h and  the  amplitude  by  F(h). 
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Figure  19.  Procedure  to  Determine  Effects  of  Nonstationary 

Excitation  on  the  Spectra  of  Atmospheric  Turbulence 
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That  part  of  the  power  spectrum  of  atmospheric  turbulence  which 
displays  a constant  -5/3  slope  on  a log-log  plot  is  referred  to  as 
the  "inertial  subrange."  As  shown  in  Appendix  B,  this  characteristic 
-5/3  slope  dictates  that  the  self-simi1ar  function  for  the  inertial 
subrange  must  be: 

F{h)  = h'/*  (66) 

This  implies  that  if  the  temporal  scale  of  a turbulence  sample  is  ex- 
panded by  a factor  of  2 and  the  amplitude  is  increased  by  a factor 
1 /3 

of  2 , the  resulting  turbulence  sample  will  have  the  same  probabilistic 

structure.  Dutton  and  Deavon  (Reference  17)  tested  for  self-similarity 
in  turbulence  samples  by  comparing  segments  of  the  same  sample  separated 
by  various  distances.  Agreement  was  very  good  for  separations  of  up  to 
50  meters.  The  failure  of  the  self-similar  concept  at  greater 
separations  can  possibly  be  attributed  to  a high  spectral  content  of 
low-frequency  waves  that  are  outside  of  the  inertial  subrange.  It  may 
be  possible  to  filter  turbulence  samples  and  then  apply  the  appropriate 
similarity  functions  for  the  various  frequency  bands. 
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The  traditional  PSD  analysis  for  calculating  the  response  of  an 
aircraft  to  continuous  turbulence  was  reviewed  with  special  attention 
to  the  underlying  assumptions.  The  degree  to  which  the  analysis  rep- 
utedly underestimates  the  number  of  occurrences  of  turbulence  in  the 
higher  load  levels  was  challenged.  It  was  shown  that  certain  miscon- 
ceptions are  responsible  for  gross  errors,  although  it  was  conceded 
that  even  a properly  applied  analysis  may  under-predict  the  number  of 
occurrences  of  extreme  turbulence.  The  error  is  reported  to  be  negligible 
for  many  applications.  Efforts  to  refine  the  procedure  by  modeling 
atmospheric  turbulence  as  a nonstationary  process  are  reviewed. 
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APPENDIX  A 

DERIVATION  OF  EXCEEDANCE  EQUATION 

This  appendix  presents  a derivation  of  the  following  equation, 
which  is  an  expression  for  the  expected  number  of  exceedances  of  a 
given  load  level  per  unit  time: 


->  -y 

N(y)  = N0  {p,eAb|  ♦ P*  tAb*} 


(67) 


where 

N = expected  number  of  exceedances  of  the  zero  load  level  per 
unit  time 

Pi  = probability  of  encountering  turbulence  of  type  i 

bi  = rms  intensity  of  turbulence  of  type  i 

A = rms  response  of  design  parameter  per  unit  rms  turbulence 
intensity 

y = load  level  of  design  parameter  for  which  exceedances  are 
calculated 


The  first  part  of  this  derivation  is  attributed  to  John  C.  Houbolt 
(Reference  11).  Turbulence  Is  assumed  to  occur  in  patches  where  each 
patch  is  a stationary  Gaussian  process.  Consider  a portion  of  the 
response  time  history  of  one  of  the  design  parameters: 


o< 

8 


Figure  20.  Typical  Response  Time  History 


Examine  a typical  crossing  of  the  load  level  yQ  in  order  to  determine 
how  much  time  is  spent  in  a narrow  band  between  yQ  and  yQ  + Ay. 
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Figure  21.  Typical  Level  Crossing 

If  the  velocity  of  the  crossing  is  y,  then  the  time  required  for  an 
upward  or  downward  crossing  is: 

(H  . fa  <«> 

In  an  arbitrary  length  of  time  , the  amount  of  time  spent  in  the  in- 
cremental area  of  hyper-space  (Ay)  (Ay)  is  by  definition: 

AT  * T,  p(y,y)  Ay  Ay  (69) 


where 


p(y,y)  = joint  probability  distribution  between  y and  y 


y * ^ * fixed  value 

The  number  of  crossings  of  the  interval  Ay  in  time  T^  will  then  be: 


= T,  M 
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The  number  of  exceedances  of  the  level  y will  be  equal  to  the  number 
of  crossings  in  one  direction  — either  with  positive  or  negative 
velocity,  but  not  both.  For  convenience,  consider  only  upward  cross- 
ings and  therefore  integrate  over  positive  values  of  velocity. 


" s T,  J y p(  y ,y ) dy 


An  aircraft  spends  only  part  of  the  time  in  turbulence  of  a particular 
type;  let  be  the  ratio  of  that  amount  of  time  to  the  total  amount 
of  flight  time: 


where 

T.j  = amount  of  time  in  turbulence  of  type  i 
T = total  amount  of  flight  ttme 

* portion  of  flight  time  in  turbulence  of  type  i 

Hence,  the  total  number  of  exceedances  of  level  y per  unit  time  is 
given  by: 


00 

"■?  i A p(vii)d'i 

* O 


Temporarily  drop  the  subscripts  i for  convenience  and  assume  that  the 
patch  of  turbulence  is  completely  Gaussian.  In  this  case  the  joint 
probability  expression  will  be: 


* Tid  * # 

2yy 
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Combine  the  two  foregoing  equations  and  consider  only  one  type  of 
turbulence: 


r y 2 Oi  .. 

N = P / — 2 « » t * dy 

i 2ir<r 


Moving  the  constants  outside  of  the  integral  sign  gives: 


p 2<r2  7 2o?” 

• * / y • * dy 

0L  CTa  «L 


2ir<ry 


In  order  to  evaluate  the  integral,  make  the  following  temporary  sub- 
stitutions: 


2 2 
y s u 


oo  * 09  4 

f.  20^  f -cu*  du 

Jy«  > dy  = / u • 

o o 


-j--  •’C‘|2  I"  * -y-  (80) 

2c  Ja  2c 
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Substituting  the  original  expression  for  c gives: 


-1 


hence 


t 2 a*  2 

Ji  • 9 d*  8 - ** 


(81) 


(82) 


Now  consider  the  power  spectrum  of  the  response  of  the  design  parameter. 
The  area  under  the  power  spectrum  is  the  mean  square  value: 


m 

O'*  * J <£(<*>)  dw 


(83) 


Assuming  that  the  process  is  random  and  stationary,  the  expected  value 
of  the  velocity  at  a particular  frequency  w is  given  by: 


y = iwy 


(84) 


Hence,  the  power  spectrum  of  the  velocity  is  given  by: 


2 

4>y[w)  * -w  <f>y (w) 


(85) 


Therefore,  the  mean  square  value  of  y in  terms  of  the  power  spectrum 
of  y is: 


i ■ -Jy*, 


(u)  d in 


(86) 


hence 


(87) 
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\ 

i 
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Substituting  this  result  Into  the  foregoing  expression  for  the  number 
of  exceedances  of  level  y gives: 


where 


(88) 


(89) 


The  above  expression  is  the  desired  result  for  flight  through  a single 
patch  of  stationary  Gaussian  turbulence  where  the  rms  value  of  the 
response  is  a . 

J 


Now  assume  that  the  aircraft  flies  through  a number  of  such  patches, 
where  each  patch  has  an  rms  response  selected  at  random  from  a 
distribution  with  a probability  density  q(oy).  The  number  of 
exceedances  of  the  level  y per  unit  time  will  be: 


N(y)  = P 


(90) 


Assume  that  q(cjy)  is  a Gaussian  distribution  that  is  restricted  to 
positive  values  of  the  random  variable,  and  therefore  the  amplitude 
has  been  appropriately  multiplied  by  2: 


«(■»)  ■ 


(91) 


where 

o = rms  value  of  design  parameter  taken  over  all  patches  of 
turbulence 


i 

j 

l 
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Combining  the  two  foregoing  equations  gives: 


(92) 


From  MATHEMATICAL  TABLES  FROM  HANDBOOK  OF  CHEMISTRY  AND  PHYSICS, 
Tenth  Edition,  page  275,  Integral  #427: 


(93) 
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> d<r,  = 


Substituting  this  result  into  the  previous  expression  for  N(y)  gives: 


N{»)  = P N0  • 


The  subscripts  i were  dropped  for  convenience  and  they  are  reinstated 


Recalling  that: 


gives 


N<»)  = N0  I P * 
i ' 


°i  s A bi 


N(y)  = H0  Z P.  • A b| 


(101) 


which  for  i = 1,  2 gives  the  desired  result. 


NCjr ) 8 N0  {p,  • Zf5i  ♦ P2  • } 
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APPENDIX  B 

SELF-SIMILARITY  RELATION  FOR  VON  KARMAN  SPECTRA 

This  appendix  demonstrates  that  the  existence  of  the  inertial 
subrange  of  atmospheric  turbulence  implies  that  the  wavelength  X and 
the  amplitude  h are  related  in  a statistical  sense  by: 


h a \ ' ' 

Consider  the  von  Karman  representation  for  the  power  spectra  of  the 
atmosphere: 


*(a)  - 


r«L  i + ■§•('•»»  lo)’ 

T [■  ♦ (1.„9La)t]"/‘ 


The  inertial  subrange  is  that  part  of  the  power  spectrum  which  displays 
a constant  -5/3  slope.  This  occurs  for  large  values  of  0: 

limit  ( 4>W\  _ ,C5/S 


limit  (Km  Q 
a-+-  \ <rz  / 


On  a log-log  plot  of  i (n)/o^  vs  S2,  the  slope  is  -5/3: 


The  power  spectral  density  function  may  be  defined  by: 


limit  3 

$(Sl)  * J h ,t  ) dt 


*(n) 
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Assume  that: 


h a ft' 


(111) 


then 


n*'-' 


■*  " Aft 


But 


Therefore 


ilSl  a ftS/S 


2y  -l  “3/3 

ft  7 oft 


hence 


X — ■ 


Substituting  into  Equation  111  gives: 

haft 


Frequency  and  wavelength  are  inversely  related.  Therefore: 


h a X' 


Equation  117  is  the  desired  result. 


(112) 


(113) 


(114) 


(115) 


(116) 


(117) 
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